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ABSTRACT 

We estimate the rate of dark matter scattering in collapsed structures throughout the 
history of the Universe. If the scattering cross-section is velocity-independent, then 
the canonical picture is correct that scatterings occur mainly at late times. The scat¬ 
tering rate peaks slightly at redshift z ~ 6, and remains significant today. Half the 
scatterings occur after z ~ 1, in structures more massive than 10^^ Mq. Within a 
factor of two, these numbers are robust to changes in the assumed astrophysics, and 
the scatterings would be captured in cosmological simulations. However, for particle 
physics models with a velocity-dependent cross-section (as for Yukawa potential in¬ 
teractions via a massive mediator), the scattering rate peaks before z ~ 20, in objects 
with mass ^ 10^ M©. These precise values are sensitive to the redshift-dependent 
mass-concentration relation and the small-scale cutoff in the matter power spectrum. 
In extreme cases, the qualitative effect of early interactions may be reminiscent of 
warm dark matter and strongly affect the subsequent growth of structure. However, 
these scatterings are being missed in existing cosmological simulations with limited 
mass resolution. 
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1 INTRODUCTION 

The cold dark matter and cosmological constant model 
(ACDM) has been successful at describing the observed 
large-scale structure. However, reported differences on 
smaller physical scales between simulations and observations 
(for a review see Weinberg et al. 2015) have raised the ex¬ 
citing question of whether one of the two main assumptions 
about the dark matter (DM), namely that it is cold with 
low thermal velocities and that it is collisionless, could need 
revising. 

Self-Interacting Dark Matter (SIDM) removes the col¬ 
lisionless assumption, such that DM particles have a cross- 
section for interacting that is sufficient to produce observ¬ 
able astrophysical effects. In the simplest model, originally 
proposed by Spergel & Steinhardt (2000), these interactions 
are elastic scatterings with an interaction cross-section that 
is independent of velocity. The DM particle collisions de¬ 
crease the central density of DM haloes and tend to make the 
DM velocity distribution isotropic, leading to more spherical 
haloes (Burkert 2000; Yoshida et al. 2000; Spergel & Stein¬ 
hardt 2000). 

Initial excitement about SIDM was related to its ability 
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to produce constant density cores in dwarf galaxies (Yoshida 
et al. 2000; Spergel & Steinhardt 2000), as well as reduce the 
amount of substructure in DM haloes. Constraints placed on 
the SIDM cross-section imply DM collisions are unlikely to 
produce significant evaporation of subhaloes, though self in¬ 
teractions can still have a noticeable effect on halo density 
profiles (Rocha et al. 2013). This could explain the results of 
a detailed comparison between local dwarf galaxies and sim¬ 
ulations made by Boylan-Kolchin et al. (2011), who found 
that the most massive DM substructures around simulated 
Milky Way-like haloes were considerably more massive than 
estimated dwarf galaxy masses made from line-of-sight ve¬ 
locity measurements (Walker et al. 2009; Wolf et al. 2010). 
If the results from the collisionless N-body simulations are 
representative of the real Universe, then there must be a sig- 
nihcant number of massive dark subhaloes around the Milky 
Way. These subhaloes that do not contain stars despite their 
large mass have been dubbed ‘too big to fail’. An alterna¬ 
tive explanation is that the most massive subhaloes do form 
stars, but that their circular velocities are below that seen 
in the collisionless DM simulations. This can be achieved 
in SIDM, where the constant density cores formed through 
DM collisions reduce the circular velocities of subhaloes. 

Since SIDM was first proposed as an alternative to col¬ 
lisionless CDM, work has been done to constrain the self- 
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interaction cross-section. Astrophysical considerations have 
included the core sizes of clusters (Yoshida et al. 2000), the 
ellipticity of clusters (Miralda-Escude 2002), evaporation of 
galaxy haloes in clusters (Gnedin & Ostriker 2001), and the 
dynamics and mass-to-light ratios of merging systems such 
as the Bullet Cluster (Markevitch et al. 2004; Randall et al. 
2008; Harvey et al. 2015). 

The tightest constraints come from galaxy cluster 
scales, where the relative velocity between DM particles is 
high. Meanwhile SIDM’s ability to solve the ‘too big to fail’ 
problem is on the dwarf galaxy scale. This was recently 
noted by Fry et al. (2015) who found that cross-sections 
consistent with cluster scale constraints could not signifi¬ 
cantly reduce the central density of haloes with peak cir¬ 
cular velocities below 30kms“^. For this reason, as well 
as the fact that many particle physics models give rise 
to them, there has been increased interest in SIDM with 
a velocity-dependent cross-section. A DM particle with a 
cross-section that decreases with increasing relative particle 
velocity (see e.g. Khoze & Ro 2014) could have an effect on 
dwarf galaxy scales where velocity dispersions are low, while 
leaving galaxy clusters relatively untouched. For this reason, 
we look at a well-motivated particle model that gives rise to 
a velocity-dependent scattering cross-section. 

Assessing the effects of dark matter particle phe¬ 
nomenology on structure formation is usually done using 
cosmological simulations. However these simulations can 
only access a hnite range of objects due to their limited reso¬ 
lution. An alternative to simulations, originally pioneered by 
Press & Schechter (1974, hereafter PS) and later extended 
by use of Excursion Set Theory (Bond et al. 1991) and con¬ 
sideration of ellipsoidal collapse (Sheth et al. 2001), is used 
to calculate ‘analytical’ mass functions. This is done using 
linear theory to evolve the density field, and assuming a sim¬ 
ple model for gravitational collapse in which regions denser 
than some density threshold collapse to form virialised ob¬ 
jects. Using the PS formalism is attractive as it allows us 
to look at all scales and redshifts simultaneously, while we 
can easily separate the contribution from haloes of different 
masses to quantities such as the mean scattering rate for 
SIDM particles through cosmic time. 

This work follows a similar procedure to Cirelli et al. 
(2009), who estimated the DM annihilation rate through 
cosmic time. The rate of interactions in a DM halo can 
be calculated given a particle model and the density pro- 
Hle of the halo. Then with a mass function (from PS theory 
or equivalent) it is possible to work out the total rate of 
scattering in the Universe. For the simplest model of parti¬ 
cle annihilation the DM cross-section, cr, multiplied by the 
relative velocity of particles, v, is constant. As the rate of 
interactions is proportional to (av) this simplifies the cal¬ 
culation relative to a case where a has some other velocity 
dependence. In this work we use DM models that have in¬ 
teraction cross-sections that differ from a oc 1/w, first using 
the simplest model for particle scattering in which cr is a 
constant. 

Our study is aimed at estimating the rate of scattering 
in DM haloes of different masses through cosmic time. The 
high redshift Universe is very dense, and were it to turn out 
that the scattering rate was therefore high, the survival of 
the first seeds of structure formation could provide a useful 
constraint on the self-interaction cross-section of DM. This 


work should also be helpful when assessing the importance 
of resolution in A^-body simulations of SIDM, because they 
can only resolve objects above a certain mass. While only 
the resolved objects from simulations are usually of interest, 
objects build up in a hierarchical fashion, such that resolved 
objects at some epoch, are made from the merging of smaller 
(potentially unresolved) objects from an earlier time. It is 
therefore important to assess whether these small objects 
that merged should have been affected by DM self interac¬ 
tions. 

The paper is organised as follows. In Section 2 we dis¬ 
cuss the calculation of the DM interaction rate through cos¬ 
mic time for a velocity-independent scattering cross-section, 
while in Section 3 we show the effects of changing the models 
and parameters that went into our original calculation. In 
Section 4 we perform the same calculation with velocity- 
dependent cross-sections, focussing in particular on two 
models recently simulated by Vogelsberger & Zavala (2013). 
Finally, we give our conclusions in Section 5. Throughout 
the paper we assume a Planck 2013 cosmology (Planck Col¬ 
laboration et al. 2014) unless stated otherwise, and also as¬ 
sume that self-interactions do not effect large scale structure 
formation. 


2 INTERACTION RATE OVER COSMIC TIME 

In this section we first discuss the number density of DM 
haloes of different masses and how this evolves with red- 
shift. By then looking at the scattering rate of DM particles 
in the haloes that exist at a particular redshift we can calcu¬ 
late the rate of DM particle scattering at that epoch. This 
calculation assumes that DM scattering is only between par¬ 
ticles within the same DM halo, and neglects the fact that 
scattering rates would be enhanced during the merging of 
DM haloes, when the relative velocities between particles 
can be larger. As haloes only spend a small fraction of time 
undergoing major mergers, the contribution of mergers to 
the integrated number of scatterings should not be too sig- 
nihcant. 


2.1 Mass function of collapsed structures 

We initially calculate the number of structures of a given 
mass using Press-Schechter (PS) theory, considering alterna¬ 
tive formulations in Section 3.2. The primordial fluctuations 
5 = {p — {p))/{p) in the Universe’s matter density field p, 
are evolved using linear theory. The spherical collapse model 
(e.g. Lacey & Cole 1993) shows that volumes of radius R in 
which the mean overdensity Sr exceeds a critical threshold 
Sr > Sc = 1.686 will collapse under their own gravity. We 
assume gravitational collapse to be immediate leading to a 
virialised halo with mass M = ^'kR^{p). 

To find these volumes, consider smoothing the density 
distribution on a scale R. Assuming the density fluctuations 
form a Gaussian random field, the fraction of the Universe 
in regions with an overdensity greater than Sc is 

This depends only on the variance of Sr on this scale. 
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Because 5r has zero mean, 


4 = ( 4 ) = D\z) I k^P{k)wUk)dk, 


( 2 ) 


where the linear growth factor D{z) governs the amplitude 
of perturbations at redshift z, and Wii{k) is the Fourier 
Transform of a real-space spherical top hat hlter of radius 
R. 

The power spectrum P{k) is obtained by multiplying 
the power spectrum of fluctuations generated by inflation 
by the Transfer Function T{k), which accounts for the dif¬ 
ferent behaviour of fluctuations that are smaller than and 
larger than the horizon during the radiation and then mat¬ 
ter dominated eras. For simplicity we use the Eisenstein & 
Hu (1998) zero-baryon CDM model in which 


Lo + Co ’ 
Lo[q) = ln(2e -I- 1.8g), 


Co(g) = 14.2 + 


731 

l + 62.5g’ 


and q is related to k by 
^ Mpc“^ 


(3) 


(4) 


where Tcmb = 2.7 ©2.7 K. We look at the effect of changes to 
the high-fc power spectrum by integrating the mass function 
down to different minimum masses, as explained in Section 

2.3. 

PS theory then interprets the fraction of the Universe’s 
volume for which 5r > 5c as the fraction of the Uni¬ 
verse’s mass that has collapsed to form objects with mass 
M ^ (p)- In this transition from smoothing over vol¬ 

umes to mass scales, it is also convenient to eliminate time 
dependence from the rms density fluctuations, i.e. we define 
the rms mass fluctuations on scale M as cm = (Jr{z) /D{z), 
such that D(z = 0) = 1. Thus the fraction of the mass in the 
Universe in collapsed objects with mass greater than M, at 
redshift z, is 


F{M,z) = 


I 


Sc/itmD{z) 




exp 



(5) 


where ^ = 5m/o-mD{z). This depends only on the rms den¬ 
sity fluctuations (in the lower limit of integration) and the 
linear growth factor. Differentiating it with respect to mass 
yields the multiplicity function 


dP 

d InM 


D) 


d In aM 


d InM 


izexp 



( 6 ) 


where we have introduced u = om'D^z) and multiplied 
by a factor of two to account for mass that is initially in 
under-dense regions/ This describes how the mass in the 
Universe is divided amongst objects of different mass and is 
plotted in the top panel of Fig. 1. 




Figure 1. Top panel. The Multiplicity Function, which shows 
how mass in the Universe is split between objects of different 
mass, as described by Press-Schechter theory. Different coloured 
lines show different redshifts. Middle panel: The interaction rate 
per particle as a function of halo mass, assuming NFW density 
profiles and the Duffy et al. (2008) concentration-mass relation, 
with a velocity-independent cross-section of lcm^g“^. Circles 
highlight the mass at which the Multiplicity Function peaks, il¬ 
lustrating a relatively constant interaction rate per unit mass in 
the Universe’s most typical haloes. Bottom panel: The product 
of curves in the two upper panels, illustrating the relative con¬ 
tribution of haloes in different logarithmic mass bins to the total 
interaction rate per particle. In this scenario, the main location 
for scatterings gradually transitions to more and more massive 
structures. 


2.2 Interaction rates in collapsed structures 

The scattering rate of an individual dark matter (DM) par¬ 
ticle 2 , with velocity Vi, is 

[ f{v')p—\vi-v'\d^v', (7) 

J rn 


^ Consider what happens when we take M ^ 0 in Equation 5. On 
small scales the rms fluctuations are very large, and the lower limit 
in the integration tends to zero. This implies F(0,z) = f, and 
only half of the mass in the Universe is in collapsed objects. On 
small enough scales the density field is always non-linear, and so 
we would expect all mass in the Universe to be in collapsed objects 


if we take M ^ 0. The missing half of the Universe corresponds to 
regions that are below the collapse threshold when smoothed on 
a scale M, but would be above the collapse threshold if smoothed 
on some larger scale. For more information, see the discussion of 
the ‘cloud-in-cloud’ problem in Bond et al. (1991). 
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where / is the velocity distribution function,^ p the local 
density, and (pI’m) the cross-section for DM-DM scattering 
(which could depend on \vi — v'\ = Upair) divided by the 
DM particle mass. Integrating over the velocity distribution 
function gives the scattering rate for a particle at position 
r, 

(8) 

For a halo of mass M containing N particles, the mean 
scattering rate per particle is 

1 ^ 

rhaic(M) =-^r,. (9) 

i = l 


Integrating over radius r gives 
1 

rhaio(M) = ^ J p{r)ri{r) dr 


1 r 

M Jo 


4^^2p2(^){^«pai^dr. 


( 10 ) 

( 11 ) 


We assume that the collapsed haloes from PS theory 
have spherically symmetric Navarro et al. (1997, hereafter 
NFW) density profiles. 


P{r) ^ ^NFW 

Pcrit (r/ra)(l + r/ra)2’ ^ ^ 

where Vs is a scale radius, 5 nfw a dimensionless character¬ 
istic density, and pcrit = SH^/S-kG is the critical density. 
We assume that the mass of a halo fills a spherical region of 
radius r 2 oo, within which the mean density is 200 pcrit and 
the total mass is M 200 . Outside this region we assume the 
density to be zero. For brevity we will also refer to r 2 oo as 
and M 200 as M. The concentration parameter is defined as 
c = r 2 oo /rs and can be related to the characteristic density 
by 


200 c® 

3 ln(l + c) — c/(l + c)' 


(13) 


Note that the NFW profile is obtained from non¬ 
interacting dark matter simulations. Dark matter scattering 
reduces the density in the centre of DM haloes, producing 
a constant density core (Burkert 2000; Yoshida et al. 2000; 
Spergel & Steinhardt 2000; Cohn et al. 2002; Rocha et al. 
2013; Zavala et al. 2013). Assuming an NFW profile, the av¬ 
erage radius at which interactions take place (assuming an 
isotropic velocity dispersion) is 0.32 independent of halo 
concentration. This is similar to the radius for which the ra¬ 
dial density profiles seen in the simulations of Rocha et al. 
(2013) first drop below the NFW prediction. These simula¬ 
tions used the maximum allowed velocity-independent cross- 
section, and so cores in other models would likely be smaller. 
Also, while the density in the centres of haloes decreases, 
DM scattering increases the velocity dispersion in halo cen¬ 
tres, which should cancel some of the effect. Nevertheless we 
acknowledge that these DM interactions are moderately self¬ 
regulating because they form cores that will tend to decrease 
the interaction rate, but proceed assuming an NFW profile 
for the DM density. If we relax this assumption then the 


^ Here / is normalised such that f f{vi)d)v\ = 1. 


scattering rates calculated would be lower, but a full treat¬ 
ment of the effect that scattering has on the phase space dis¬ 
tribution of haloes, and so the subsequent scattering rates, 
requires full N-body simulations that are beyond the scope 
of this paper. 

To calculate the mean pairwise velocity of particles, we 
integrate over the velocity distribution functions of particle 
pairs. Assuming that their velocities are isotropic and follow 
a Maxwell-Boltzmann distribution^ with one dimensional 
velocity dispersion uid, this gives (upair) = (4/v^)criD. For 
an NFW halo, the velocity dispersion of particles is (Lokas 
& Mamon 2001) 


TT^ — In(cs)- — 

cs 

_ 1 _ 6 _ A _^ _ 2 \ 

(1 -I- cs)^ 1 -|- cs \ c^s^ cs 1 -I- cs/ 

X ln(l-I-cs)-I-31n^(l-I-cs)-I- 6 Li 2 (—cs) , 

(14) 

where s = r/r 2 oo, g{c) = [ln(l -h c) — c/(l -h c)]“^, and 
Li 2 (a;) is the dilogarithm (commonly referred to as Spence’s 
function), defined by 


2 / N . -.2 GM 

o-id(s,c) = -c 5 r(c)s(l-h cs) - 

2 Tv 


Li2(a;) 



ln(l — u) 
u 


dw. 


(15) 


Returning to Equation (11), and changing integration 
variable from r to s, we find 


rhaio(M,rv,c) = 16v^^— / s^p^(s, c) o-id(s, c) ds, (16) 

M m Jo 

where we have now assumed that the DM-DM cross-section 
is velocity-independent (this restriction is relaxed in Sec¬ 
tion 4). Both p(s, c) and ctid(s, c) depend on the virial mass 
and radius of a halo, and can be written as dimensionless 
functions of s and c multiplied by the dimensional quanti¬ 
ties M/r® and -y/ GM/r^ respectively. We can then see that 
Fhaio will be a function of the halo concentration scaled by 
power-laws in M and ry. Specifically, at fixed cross-section 
and halo concentration, Fhaio oc 

At a particular cosmic time, M = M 200 and = ^200 
are not independent, because M 2 oo/r|oo oc pcrit(z) by defi¬ 
nition. Using this, we find Fhaio oc pA®, such that 


rhaio(M, 2 , (cr/m),c) = rhaio(Mo,2o, (cr/m)o,c) 


/ M 


^ y Pcrit (g) y (cr/m) \ 

VPcrit(zo)y \{o/m)o)' 

(17) 

We calculate rhaio(Afo, 20 , (u/m)o, c) with Mo = IO’^’^Mq, 
2 o = 0 and (aj’mjo = 1 cm^ g“^, by numerically integrating 
Equation (16). We can then calculate Fhaio for haloes with 
different masses and at different redshifts using Equation 
(17). 


^ High resolution simulations of CDM report departures from 
Gaussianity for the distribution of velocity components along the 
principal axes of the velocity dispersion tensor (Vogelsberger et ah 
2009), but this approximation is sufficient for our work. 
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At fixed mass, redshift and cross-section, Fhaio is found 
to increase significantly with increasing halo concentration. 
The logarithmic slope of the rhaio(c) relation is ~ 1.7 for 
c = 5, and ~ 2.5 for c = 30, with Thaio oc (? for con¬ 
centrations around 10. As halo concentrations generally de¬ 
crease with increasing halo mass, the mass dependence of 
Thaio is suppressed below the Thaio ric seen in Equa¬ 

tion (17). The overall form of Thaio(M, z) depends upon the 
concentration-mass-redshift relation. Following Duffy et al. 
(2008, hereafter DOS), we shall initially assume 

/ M \ -0 081 

Using this c{M, z) relation we show Thaio (AT, z) in the middle 
panel of Fig. 1. Thaio increases rapidly with increasing red- 
shift at fixed mass, and increases with mass at fixed redshift. 
As objects grow in mass through cosmic time, the scattering 
rate in typical haloes at each redshift evolves slowly. Note 
that several more recent works show that this simple power 
law dependence of c{M) should flatten at low masses, as is 
discussed in Section 3.1. 

2.3 DM’s cosmic scattering rate 

Multiplying the multiplicity function from Section 2.1 by 
the interaction rate in individual haloes from Section 2.2 
gives the contribution of haloes of different mass to the total 
rate of particle scattering in the Universe (see bottom panel 
of Fig. 1). Integrating this quantity over all halo masses 
at different redshifts yields the mean scattering rate of all 
particles at that redshift, T(«), which we refer to as the 
‘Cosmic Scattering Rate’. This is plotted in Fig. 2, where it 
can be seen that after a gradual rise from the early universe 
to z « 6, T( 2 ;) is constant to within a factor of two to the 
present day. 

For this analysis, we assume that haloes form down to 
masses of 10“^^ Mq. In the real Universe, self-interacting 
dark matter creates a small-scale cut-off in the power spec¬ 
trum due to collisional damping. For DM composed of 
weakly interacting massive particles (e.g. neutralinos), the 
minimum mass of collapsed objects is ~ 10“® Mq (Hofmann 
et al. 2001). If DM were axions then this minimum mass 
would be ~ Mq (Kolb & Tkachev 1996). For the gen¬ 

eral class of self-interacting dark matter models that lead to 
astrophysically interesting scattering rates in the late-time 
Universe, collisions in the early Universe suppress power on 
larger scales, or even introduce acoustic oscillations in the 
dark matter-dark radiation system (Buckley et al. 2014). 
There is a rich possible phenomenology affecting the high-fc 
power spectrum, which is highly model-dependent. 

We investigate the approximate effect of a cutoff in the 
power spectrum by integrating Thaio(dF’/dln AT) down to 
different minimum masses, ATmin, shown as the extra lines in 
Fig. 2. Furthermore, in numerical simulations, only haloes 
above a given mass scale are resolved, and only the DM 
interactions above those scales can be tracked. We therefore 
include lines with large ATmin in Fig. 2, to act as predictions 
for the expected scattering rate in cosmological simulations. 
Note that the results as ATmin —> 0 converge particularly 
slowly for the DOS concentration-mass relation, due to the 
high concentration of very small haloes. Nevertheless, these 
results are less sensitive to changing ATmin than those for a 



Figure 2. Top panel: The mean scattering rate of particles in 
the Universe calculated from Press-Schechter theory, assuming 
the NFW density profile, the DOS concentration-mass relation 
and cr/m = fcm^g”^. The different lines count only scatterings 
in haloes more massive than IO^^'Mq (bottom line) to 10 “^^Mq 
( top line). The scattering rate varies by less than a factor of two 
from 2 Ri 6 onwards. Bottom panel. The mean cumulative num¬ 
ber of interactions that particles have undergone as a function 
of redshift. The different lines again include only those inter¬ 
actions in haloes more massive than a given threshold. With a 
velocity-independent cross-section, most scattering is at late red- 
shifts where there is more time. This results in most scattering 
being in high-mass haloes, so that Nscatteriz = 0) varies by less 
than 25% between Mmin = 10 “^^Mq and ATmin = IO^'^Mq. 


simple annihilation channel where crUpair is constant (Mack 
2014) and low mass haloes make a dominant contribution to 
the total scattering rate. 

In addition to the rate of cosmic scattering, an interest¬ 
ing quantity is the mean cumulative number of interactions 
that particles have undergone. As each scattering event is 
a two-body interaction, this is twice the number of inter¬ 
actions per particle. We call this quantity Nscatter and plot 
it as a function of redshift in the bottom panel of Fig. 2. 
While the cosmic scattering rate is markedly different at 
intermediate and high redshifts when using different min¬ 
imum masses, the values of Afat:atter(z = 0) are more ro¬ 
bust. For (n/m) = Icm^g”^, A'scatter(2 = 0) is 0.87 with 
ATmin = 10“^^ Mq and 0.68 with ATmin = 10^° Mq. 

A noticeable feature of T( 2 ) in the upper panel of Fig. 
2 is the upturn after z ~ 1. This is not present when us¬ 
ing more recent c(AT, z) relations with more complex red¬ 
shift dependences than the simple (1 -b 2 )“°'^'^ in the DOS 
relation. This upturn is not physical, and arises because 
the concentration is defined in terms of r 2 oo which in turn 
depends on pcrit. When the Universe is matter-dominated 
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Pcrit oc ( 1 + 2 )^, such that at fixed halo mass r 2 oo oc (1 + 2 )“^. 
At late times, when there is a significant dark energy contri¬ 
bution to the Universe, the evolution of pcrit slows and is no 
longer given by a simple power law in (1 -|- 2 ). This affects 
the r 2 oo of haloes, and hence halo concentrations, such that 
a simple power-law cannot accurately capture c(M, 2 ). 

3 SENSITIVITY TO ASTROPHYSICAL 

ASSUMPTIONS 

In the previous section we considered the redshift depen¬ 
dence of DM scattering rates and showed that with a 
velocity-independent cross-section, the mean rate of particle 
scattering in the Universe initially grows and then starts to 
decrease after 2 ~ 6, dropping by less than a factor of two 
to the present day. In this section we explore the sensitivity 
of this result to the assumptions of the model. 

3.1 Concentration-mass-redshift relations 

The concentration-mass-redshift relation, c(M, 2 ), of DOS is 
attractive for its simplicity and because over a small range 
of redshifts and halo masses, concentrations can be well fit 
by simple power laws in M and (1 -|- 2 ). However, numeri¬ 
cal studies that have resolved structures over a wide range 
of halo masses have found that concentrations are not well 
fitted by simple power laws. Examining the results of the 
Millennium Simulation (Springel et al. 2005) from 2 = 3 
to 2 = 0 it is clear that the form of c(M, 2 ) is not sep¬ 
arable, with the mass dependence evolving with redshift 
(Gao et al. 2008). This evolution takes the form of a flat¬ 
tening of the c(M) relation at increasing redshift, such that 
concentrations of very massive galaxy cluster haloes evolve 
only weakly with redshift while the concentrations of smaller 
haloes decrease rapidly with increasing redshift. 

The c(M, 2 ) relation is found to be remarkably complex, 
particularly when considering the dependence on cosmolog¬ 
ical parameters. Prada et al. (2012, hereafter P12) find that 
this complex relationship is a result of the ‘wrong’ physical 
quantities, M and 2 , being used. Analogous to studies of the 
halo mass function, in which a much simpler fitting formula 
is possible when one considers the mass function as a func¬ 
tion of Incr)^/ rather than a function of M (Jenkins et al. 
2001), the c(ln(jM) relationship is found to be simpler than 
c{M). 

The behaviour of this relationship can be explained by 
models in which the concentration of a halo is related to its 
accretion history (Wechsler et al. 2002; Zhao et al. 2003). 
Ludlow et al. (2014, hereafter L14) found that if the mass 
of a halo, M{z), was plotted against the critical density, 
Pcrit ( 2 :), then the relationship M(pcrit) was well fit by an 
NEW profile, with associated concentration cmah- They also 
found a simple relation between cmah and the concentra¬ 
tion of the halo, allowing the concentration-mass relation to 
be predicted from the mass-accretion history of haloes. The 
statistics of the mass-accretion of DM haloes can be found 
from simulations, or calculated using the conditional proba¬ 
bilities^ found in extensions of PS theory (Bond et al. 1991; 
Bower 1991; Lacey & Cole 1993; Kauffmann & White 1993). 

^ The conditional probability that the material making up an 



Figure 3. The cosmic scattering rate calculated using the 
concentration-mass-redshift relations of Duffy et al. (2008), Prada 
et al. (2012), Ludlow et al. (2014), Dutton & Maccio (2014), 
Diemer & Kravtsov (2015), and Correa et al. (2015). These were 
calculated assuming a Planck 2013 cosmology (Planck Collabora¬ 
tion et al. 2014), a PS mass-function, and tr/m = lcm^g“^, 
counting the contribution from all haloes more massive than 
10“^^ Mq. Lines become dashed for redshifts where authors state 
their relationships may not be valid. 

Different methods for measuring c(M, 2 ), either from 
simulations or analytical calculations, give similar results 
around the peak of the multiplicity function (M « M*), 
but differ significantly at high and low masses. While the 
cosmic scattering rate is dominated by haloes around M*{z), 
the scattering rate in haloes is highly sensitive to the halo 
concentration, and so even small differences between c(M, 2 ) 
relations can lead to significant changes in r( 2 ). In Fig. 3 we 
show r( 2 ) calculated as in Fig. 2 but for a variety of c(M, 2 ) 
relations. 

Noticeable in Fig. 3 is that using c(M, 2 ) from L14 gives 
a scattering rate at intermediate redshifts a factor of two 
above that found using other c(M, 2 ) relations. The L14 
analytical model was calculated for relaxed haloes, which 
are generally dynamically older, making them more concen¬ 
trated than unrelaxed haloes of a similar mass. The cuts 
made to remove unrelaxed haloes are one of the two main 
reasons why c{M, 2 ) relations from simulations disagree with 
each other, the other being the way in which c is calculated 
from a mass distribution. For example, Prada et al. (2012) 
calculate c from the ratio Umax/V 2 oo, where Umax and U 200 
are the maximum circular velocity and the circular velocity 
at r 2 oo respectively, while Diemer & Kravtsov (2015) find c 
by directly fitting the radial density with an NEW profile. 

3.2 Mass Function Prescription 

It is known that the PS formula does not provide an exact 
fit to the mass function from simulations. Specifically, it un¬ 
derestimates the number of rare objects in the ‘high-mass 

object of mass Mi at redshift 2 i is in an object of mass Mq at 
redshift zq. 
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tail’, with an overestimate of the amount of mass around 
the peak of the multiplicity function (see e.g. Jenkins et al. 
2001). A better fit to the mass function from simulations 
was achieved by Sheth & Tormen (1999, hereafter ST), who 
found that compared to PS, Equation (6) becomes: 

[l + i^exp - 

(19) 

with A = 0.3222, a = 0.707 and p = 0.3. We note that our 
definition of v is different from that in the ST paper, with 
usT = The original PS formula can also be described by 
Equation (19) with A = 0.5, a = 1 and p = 0. 

The ST mass function increases the number density of 
the most massive objects compared to the PS mass function, 
providing a better fit to simulations (see e.g. Reed et al. 
2007). While these differences can be extremely important 
for some studies (e.g. counting the number density of mas¬ 
sive clusters) we find that the different mass functions do 
not have a large effect on our results. This is because the 
scattering rate per unit mass in DM haloes increases only 
gently with increasing halo mass, as can be seen in the mid¬ 
dle panel of Fig. 1. The shape of r(2) is similar when either 
a PS or ST mass function is used, although the normalisa¬ 
tion is slightly lower for the latter. By redshift zero there are 
~ 20% fewer DM interactions with an ST mass function. 


dF 

d InM 


= A, 


d In aM 
dlnM 


3.3 Varying Cosmological Parameters 

Similar to changing the formalism used to calculate the 
multiplicity function, small changes to the Cosmological 
Parameters leave the cosmic scattering rate relatively un¬ 
changed because of the weak mass dependence of rhaio(Af). 
Using c(M, z) from DOS, we find that changing cosmologi¬ 
cal parameters from Planck 2013 to WMAP9 decreases the 
mean number of interactions per particle by redshift zero, 
Nscattei{z = 0), by 12%. This is driven by Planck’s larger 
value for Dm, resulting in larger critical densities at early 
times. Using earlier WMAP results leads to similar changes, 
except for WMAP3 for which the anomalously low Dm and 
erg lead to a 33% reduction in Ascatter(z = 0). 

The concentration-mass-redshift relation also depends 
on cosmological parameters, which is made explicitly clear 
by relations that relate c to ctm rather than M directly 
(e.g. Prada et al. 2012; Diemer & Kravtsov 2015). This cos¬ 
mology dependence of c{M, z) makes little difference when 
moving from Planck 2013 to WMAP9, but further reduces 
the scattering rate for a WMAP3 cosmology such that 
Vacatter(2 = 0) is 40% lower than with a Planck 2013 cos¬ 
mology, using c{M,z) from P12. This increased difference, 
beyond that seen for a cosmology independent c{M, z), can 
be understood by noting that haloes of a particular mass 
form later with smaller as, and are therefore less concen¬ 
trated. 


3.4 Scatter in the Concentration-Mass Relation 

So far we have assumed that given the mass of a halo 
we know its concentration through the concentration- 
mass relation. In practice this relation has some scat¬ 
ter around it, which will impact on the mean scattering 


rate of haloes. From Equation (17) the concentration de¬ 
pendence of the scattering rate in haloes is described by 
rhaio(Mo, 2 o, ((7/m)o, c). TMs is non-linear in c, such that 
even symmetric scatter in c at fixed mass will alter the mean 
scattering rate in haloes of that mass. 

To discuss how Fhaio is affected by scatter in c, it will be 
useful to introduce cq, the value of c implied by the c(M, z) 
relation. Dolag et al. (2004) find that for haloes of fixed mass 
and redshift, Inc is normally distributed. If we assume that 
Inc is normally distributed with mean Inco and variance 
af^^, then c follows a log-normal distribution, with proba¬ 
bility density function 


P(c) = 


CCTlr 




exp 


(Inc — Inco)^ 
2 a? 

In c 


( 20 ) 


Including a log-normal distribution of concentrations 
at fixed mass and redshift leads to an increase in P at all 
concentrations, related to the long tail of the distribution 
towards high values, as well as a shift in the expectation 
value of c.® If c(M, z) in DOS was a measure of the mean c 
for a particular mass of halo, then we would have to make 
the change Inco —>■ Inco — (jf„,,/2 in Equation (20) to keep 
(c) = Cq. However, the c(M, z) relation in DOS was acquired 
by fitting to the median values of c in each mass bin at each 
redshift. The median value of c from the probability density 
function in Equation (20) is simply exp(lnco) = co as re¬ 
quired. Dolag et al. (2004) found that nine ~ 0.22, almost 
independent of the cosmological model. This corresponds to 
a standard deviation in logjQ c of 0.1, or a scatter of 0.1 dex. 
We find that the shape of r( 2 ) is effectively unchanged by 
scatter in c(M, z), but that the normalisation increases with 
increasing scatter. For a 0.1 dex scatter, the normalisation 
increases above that of the scatter-free case by less than 
15%. 


4 VELOCITY-DEPENDENT 
CROSS-SECTIONS 

Having calculated the rate of DM scattering through 
cosmic time assuming that the cross-section is velocity- 
independent, we now lift this assumption, and perform the 
same calculation with velocity-dependent DM-DM cross- 
sections. 


4.1 Particle model 

For velocity-dependent cross-sections we use the vdSIDMa 
and vdSIDMb models described in Vogelsberger & Zavala 
(2013). These are well-motivated by particle physics, and 
describe the transfer cross-section for elastic scattering me¬ 
diated by a new gauge boson of mass m^. This results in 
an attractive Yukawa potential with coupling strength Qc. 
These interactions are analogous to the screened Coulomb 
scattering in a plasma, for which the momentum-transfer 


® For the distribution in Equation (20), the expectation value of 
c is given by (c) = exp (incp + cr'^^/2) > eg. 
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cross-section can be approximated by 


4.2 vdSIDM cosmic scattering rates 


(JT 


477 

22.7 

13^ ln(l + /3-i). 

P < 0.1 

877 

22.7 

13^ (1 -hi.5/31-®®) 

0.1 < /3 < 10® 

77 

22?7 


, ,3 >10®, 



(21) 

= 

TTUmax/Vpair and 

= 22.7lml (Feng 


et al. 2010; Finkbeiner et al. 2011; Loeb & Weiner 2011). 
Here Vmax is the velocity at which (or Vpair) peaks, with 
o'T(vms.x) ~ We have also introduced the “momentnm- 
transfer cross-section”, ctt, defined as 


ar = J{1 — cos6)'^{0) dH 


( 22 ) 


/ I 1 

(1 — cos 0)^(0) d cos 0, (23) 

-1 <Ail 

where is the differential cross-section, assumed to be az- 
imuthally symmetric, which describes the probability of par¬ 
ticles scattering into a patch of solid angle dH. The trans¬ 
fer cross-section is an effective scattering cross-section that 
is useful in describing angularly dependent cross-sections 
(where ^ is not constant). For isotropic scattering (^ = 
constant) the transfer cross-section is simply ax = cr, while 
in general the mean momentum transfer for a scattering 
process with transfer cross-section or is equal to the mean 
momentum transfer for isotropic scattering with a — ax- 
Throughout the rest of this paper, when calculating the rate 
and number of particle scattering events we will use ox as 
if it were the cross-section i.e we will calculate an effective 
rate of particle scatterings that is the rate of isotropic scat¬ 
tering events that would lead to the same rate of momentum 
transfer.® 

The velocity-dependent cross-section in Equation (21) 
leads to noticeable changes in rhaio(Af). The cross-section 
diverges as the pairwise velocity tends to zero, such that 
scattering in low mass haloes (with typical velocities less 
than Umax) is enhanced above the constant cross-section 
case. For Upair > Umax, Ox oc leading to a strong sup¬ 
pression of the scattering rate in DM haloes with velocity 
dispersions larger than Umax- 

The vdSIDMa and vdSIDMb models have values of 
Ox’^^l'm and Umax chosen to maximise the self-interaction 
rate at the typical velocity dispersion of Milky Way dwarf 
spheroidals, while avoiding known astrophysical constraints 
on the cross-section. Specifically, vdSIDMa and vdSIDMb 
have Umax = 30kms“^ and = 3.5cm^g“^, and 

Umax = 10kms“^ and o^^^jm = 35cm^g“^ respectively. 


® In general, particle orbits within a DM halo are approximately 
isotropic, so there is no preferred direction for particle scattering. 
In these cases, the momentum transfer cross-section accurately 
captures the effects of scattering. However, this may not be the 
case for systems where there is a preferred direction along which 
particles approach (Kahlhoefer et al. 2014), such as in a merger. 


The calculation of the DM scattering rate F(z) proceeds in 
a similar manner to Section 2, in that we first find the dis¬ 
tribution of haloes of different mass ( ^ then find 

the scattering rate per unit mass in these haloes, rhaio(M). 
However, the calculation of rhaio(A^) is complicated by the 
velocity-dependent cross-section, because the cross-section 
can no longer be taken outside the integral in Equation 
(11). Instead, we find (o'Upair)(u) by numerically integrating 
o'T(upair) Upair Over the probability distribution of pairwise 
velocities, again assuming that the velocities of individual 
particles are drawn from a Maxwell-Boltzmann distribution 
function with ID velocity dispersion ctid. This yields 

(fjTu) (uid) = 3 ^ _ [ (TT(u)u®e“’' du. (24) 

2crfDV^ J 

Then with o'id(u) from Equation (14) we can find (oxv) (r), 
which we use in the numerical evaluation of Equation (11) 
to calculate Fhaio(M). Combining Fhaio(A4) with the multi¬ 
plicity function we can calculate the cosmic scattering rate 
as in Section 2.3. 

The scattering rate through cosmic time is plotted 
for vdSIDMa and vdSIDMb in Fig. 4. In contrast to the 
velocity-independent case in Fig. 2, the scattering rate is 
now displayed on a logarithmic scale. It peaks at redshift 
20 — 30 and falls by two orders of magnitude before 2 = 0. 
Most interactions thus occur at early times as can be seen in 
Fig. 5. Half occur before 2 = 5.7 for vdSIDMa and 2 = 7.2 
for vdSIDMb (in the latter case, the Universe is ~ 5% of its 
present age). This is in stark contrast to the gentler evolu¬ 
tion of r( 2 ) with a constant cross-section (c.f. Fig. 2), where 
half the interactions occur after 2 = 0.96. 

With a velocity-dependent cross-section, most scatter¬ 
ings also occur in low mass haloes with typical velocities 
V < Umax- Raising the minimum mass of considered haloes 
Mmin from 10“^^ Mq to 10® M 0 lowers the number of inter¬ 
actions by redshift zero by a factor of six, which can be seen 
in Fig. 6 (for which we introduce No = A^scatter( 2 : = 0)). For 
the constant cross-section case, the same change leads to a 
decrease in AIscatter( 2 : = 0) of only 10%. 

The choice of concentration-mass-redshift relation be¬ 
comes more important when the cross-section is velocity- 
dependent, because different c(M, 2 ) disagree most for low 
mass haloes and at high redshift. In particular, the sim¬ 
ple power law relation from DOS predicts low-mass haloes 
to be much more concentrated than more recent relations 
in which c{M) flattens at low mass. This recovers (a less 
extreme version of) what is seen in estimates of the DM 
annihilation rate, where (oVpair) is usually assumed to be 
constant, resulting in an even larger fraction of interactions 
occurring in low mass haloes, and hence a cosmic scattering 
rate with strong dependence on c(M, 2 ) (Mack 2014; Correa 
et al. 2015). 

As well as the three particle models already discussed 
(velocity-independent, vdSIDMa and vdSIDMb), we include 
in Fig.s 5 and 6 plausible but more extreme velocity- 
dependent models with lower Umax. We need not specify 
the normalisation of Ox’^^/m for these calculations, but it 
can be chosen to solve small-scale problems at dwarf galaxy 
scales, while eluding constraints at cluster scales. As Umax 
is lowered, a larger fraction of interactions happen at high 
redshift and in low-mass haloes. For the most extreme case 
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vdSIDMa vdSIDMb 



Figure 4. Scattering rates (top row) and cumulative number of scatters (bottom row) as a function of redshift, for two different 
velocity-dependent scattering cross-sections. The left column is for vdSIDMa which has Umax = 30kms“^ and fTmax/m = 3.5cm^g“^; 
while vdSIDMb (right column) has Umax = 10kms“^ and (Tmax/m = 35cm^ The different line colours correspond to different values 
for Mmin of 10®, 10^, 1,10“^, 10“®, and 10“^^ Mq, with both F and Aigcatter monotonically increasing with decreasing Mmin- The solid 
lines are for the DOS concentration-mass-redshift relation, while the dashed lines use the P12 c{M,z). Unlike the constant cross-section 
case in Fig. 2, r( 2 ) is now plotted on a logarithmic scale as the scattering rate is larger by around two orders of magnitude at high 
redshift compared to redshift zero. 


considered, with Umax = 10“® km , half of the interactions 
have occurred by 2 = 19, and half occur in haloes of mass 
< 10“® Mq. We stress that such models cannot be excluded 
on particle physics grounds, but it is unclear whether the 
large number of scatterings in such low mass haloes would 
leave a detectable signal in the present day universe. 


5 CONCLUSIONS 

We have presented an analytical calculation of the mean 
rate of dark matter-dark matter scattering events, for parti¬ 
cle physics models with a velocity-independent or velocity- 
dependent cross-section. In all our calculations, we assume 
that the self-interactions are a small perturbation to ACDM 
and do not, for example, change the overall growth of struc¬ 
ture. 

For particle physics models with a velocity-independent 
interaction cross-section, our results match the canonical 
picture in which most scatterings occur in massive struc¬ 
tures Sb 10^^ Mq at late times z ^ 1. Our calculations are 
found to be robust to current uncertainties in cosmological 
parameters as well as variations in the mass function used. 
They are also insensitive to the high-fe power spectrum (be¬ 
cause most scattering events occur in haloes more massive 
than the cut-off scales due to DM self-interactions in the 


early Universe). The main source of uncertainty in the re¬ 
sults is the concentration-mass-redshift relation c{M,z). Its 
unknown form at high redshift and low mass propagates into 
a factor of almost three discrepancy in the scattering rate 
at intermediate redshifts (2 « 10). However, the scattering 
rate changes by only a factor of two over most of cosmic 
time, and different c{M, 2 ) relations give similar results af¬ 
ter 2 « 1, where there is more time. Consequently, the total 
number of interactions during the entire history of the Uni¬ 
verse is uncertain to only a factor of ~ 2. 

For particle physics models with a well-motivated ve¬ 
locity dependence, the scattering takes place mainly in low 
mass objects 10 ^ Mq at early times 2 7. The scatter¬ 

ing rate r( 2 ) peaks at earlier redshifts 2 ~ 20, with a pro¬ 
nounced peak two orders of magnitude higher than the scat¬ 
tering rate at the present day. These numbers are more sen¬ 
sitive to the choice of cosmological and astrophysical pa¬ 
rameters, and are dominated by regimes in which the mass 
function and concentration-mass-redshift relation are least 
well known. This minimum mass of considered structures, 
Mmin, is particularly important, with changes to the small- 
scale power spectrum induced by DM scattering affecting 
the cosmic scattering rate. 

The dominance of high redshift scatterings in velocity- 
dependent models profoundly changes their influence on the 
evolution of structure, and may alter the best strategy to 
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Figure 5. When do scatterings happen? The cumulative number 
of interactions as a function of redshift, normalised to unity at 
redshift zero. The different colours correspond to different parti¬ 
cle models for the DM, while the solid and dashed lines are for 
the DOS and P12 c{M,z) relations respectively. All curves were 
calculated using Mmin = 10“^^ Mq. The number in brackets in 
the legend is Ascatter(-2 = 0) for the relevant model. These are not 
present for the models with specified Umax, which represent vd- 
SIDM models with unspecified Velocity-dependent models 

with low Umax lead to more interactions in haloes with low in¬ 
ternal velocities, pushing scattering towards high redshifts where 
collapsed objects are less massive. 

search for observational signatures. DM particle interactions 
lead to a transport of particles away from the centres of 
structures (Kochanek &; White 2000), replacing the cusps 
found in collisionless CDM simulations with constant den¬ 
sity cores.' If the SIDM interactions are effectively confined 
to high redshift, then they may lead to a smearing of small- 
scale structure more qualitatively reminiscent of warm dark 
matter. The affected DM structures are also the hosts of the 
first galaxies, and it is interesting to consider what impact 
cored haloes could have on early galaxy formation. 

High redshift scattering in low-mass objects also has 
important consequences for attempts to simulate vdSIDM 
cosmologies. Most scatterings occur in low mass haloes at 
high redshift that would not be resolved in typical cosmo¬ 
logical simulations, but the unresolved interactions could 
be important for the later dynamics of particles. The large 
number of self-interactions would create DM cores in high- 
redshift haloes, and it then becomes an important question 
- on which there seems little consensus - whether or not 


^ Haloes in which there have been a larger number of interac¬ 
tions presumably have larger cores. However, we caution against 
qualitative attempts to determine the scattering rate from core 
sizes: estimates of the core size and ellipticity of a galaxy cluster 
halo (Miralda-Escude 2002) overestimated the effect of SIDM by 
a factor 50 compared to full simulations (Peter et al. 2013). Fur¬ 
thermore, strong gravitational lensing measurements of the cores 
in low redshift clusters (Sand et al. 2003; Newman et al. 2013) are 
subject to projection effects. Particularly when SIDM lowers the 
central density, material at large radii significantly contributed to 
the 2D projected density. 



Figure 6. Where do scatterings happen? The fraction of scat¬ 
terings by redshift zero that occur in haloes more massive than 
A/min, normalised to unity for Mmin = 10“^^ Mq. Different line 
styles are as in Fig. 5, with colours corresponding to particle 
models, and solid or dashed lines representing the DOS or P12 
concentration-mass-redshift relations respectively. Models with 
velocity-independent cross-sections have more of their scatterings 
in high-mass haloes compared to velocity-dependent cases, where 
the typical halo mass in which most interactions happen is an 
increasing function of Umax- 


the mergers of small cored haloes form cores that persist in 
large haloes at the present day. 
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